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Abstract 


The upwind leapfrog or Linear Bicharacteristic Scheme (LBS) has previously been implemented and 
demonstrated on electromagnetic wave propagation problems. This report extends the Linear Bicharac- 
teristic Scheme for computational electromagnetics to model lossy dielectric and magnetic materials and 
perfect electrical conductors. This is accomplished by proper implementation of the LBS for homogeneous 
lossy dielectric and magnetic media and for perfect electrical conductors. Heterogeneous media are modeled 
through implementation of surface boundary conditions and no special extrapolations or interpolations at 
dielectric material boundaries are required. Results are presented for one-dimensional model problems on 
both uniform and nonuniform grids, and the FDTD algorithm is chosen as a convenient reference algorithm 
for comparison. The results demonstrate that the explicit LBS is a dissipation-free, second-order accurate 
algorithm which uses a smaller stencil than the FDTD algorithm, yet it has approximately one-third the 
phase velocity error. The LBS is also more accurate on nonuniform grids. 


1 Introduction 

Numerical solutions of the Euler equations in Computational Fluid Dynamics (CFD) have il- 
lustrated the importance of treating a hyperbolic system of partial differential equations with the 
theory of characteristics and in an upwind manner (as opposed to symmetrically in space). These 
two features provide the motivation to use the Linear Bicharacteristic Scheme (LBS), or the up- 
wind leapfrog method, for the construction of many practical wave propagation algorithms. In a 
hyperbolic system, the solutions (i.e. waves) propagate in preferred directions called characteris- 
tics. A characteristic can be defined as a propagation path along which a physical disturbance is 
propagated [1]. The relevance to Maxwell's equations is intuitively obvious because electromag- 
netic waves have preferred directions of propgation and finite propagation speeds. 

The upwind leapfrog method has a more compact stencil compared with a classical leapfrog 
method. Clustering the stencil around the characteristic enables high accuracy to be achieved 
with a low operation count in a fully discrete way [2]. This leads to a more natural treatment of 
outer boundaries and material boundaries. The LBS treats the outer boundary condition naturally 
without nonreflecting approximations or matched layers. The interior point algorithm predicts 
the outgoing characteristic variables at the domain boundaries. Through knowledge of the wave 
propagation angle, the local coordinates can be rotated to align with the characteristics, at which 
the boundary condition becomes almost exact. Therefore, no extraneous boundary condition or 
matched layers are required, which can introduce errors into the solution. The LBS also offers a 
natural treatment of dielectric interfaces, without any extrapolation or interpolation of fields or 
material properties near material discontinuities. 

The LBS was originally developed to improve unsteady solutions in computational acoustics 
and aeroacoustics [3]-[8]. It is a classical leapfrog algorithm, but is combined with upwind bias in 
the spatial derivatives. This approach preserves the time-reversibility of the leapfrog algorithm, 
which results in no dissipation, and it permits more flexibility by the ability to adopt a charac- 
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teristic based method. The use of characteristic variables allows the LBS to treat the outer com- 
putational boundaries naturally using the exact compatibility equations. The LBS offers a central 
storage approach with lower dispersion than the Yee algorithm, plus it generalizes much easier to 
nonuniform grids. It has previously been applied to two and three-dimensional free-space elec- 
tromagnetic propagation and scattering problems [4], [7], [8]. 

The objective of this report is to extend the LBS to model both homogeneous and heteroge- 
neous lossy dielectric and magnetic materials and perfect electrical conductors (PECs). Results 
are presented for several one-dimensional model problems, and the FDTD algorithm is chosen 
as a convenient reference for comparison. Sections 2 and 3 present the LBS implementation for 
homogeneous and heterogeneous materials, respectively. Section 4 discusses the outer radiation 
boundary condition and Section 5 reviews the Fourier analysis. Finally, Section 6 presents results 
for one-dimensional model problems and Section 7 provides concluding remarks. 


2 Homogeneous Materials 


Maxwell's equations for linear, homogeneous and lossy media in the one-dimensional TE case 
(taking d/dy = djdz = 0) are 
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where a and a* are the electric and magnetic conductivities, respectively. Using the electric dis- 
placement D = eE and making the substitution c = 1/077 gives 
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The procedure for the LBS is to transform the dependent variables D y and H, to characteristic vari- 
ables. Based upon the numerical method of characteristics for electromagnetics [9], information 
propagates (in one space dimension) along the characteristic curves specified by the characteristic 
equations 
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which show that waves propagate in the ±x directions with a finite physical speed of c. The 
algorithm developed here is the simplest leapfrog scheme described by Iserles [10] combined with 
upwind bias, or simply, the Linear Bicharacteristic Scheme (LBS). To transform (3) and (4) into 
characteristic form, we first multiply (4) by c and then add and subtract from (3) to give 
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Now define 

P = Dy + -H z (8) 

c 

Q = Dy--H z (9) 

C 

to represent the right and left propagating solutions, respectively. P and Q are otherwise known 
as the characteristic variables. Using these definitions, (6) and (7) can be rewritten as 
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It is convenient to define and store the collowing coefficients before 
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Equations (10) and (11) can be rewritten more simply as 
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To develop the discretized algorithm for a one-dimensional system, the stencils of Figure 1 are 
proposed for the LBS. To solve the wave propagation problem without introducing dissipation, it 
is necessary that the stencil have central symmetry so the scheme employed is reversible in time 
[2]. The stencil in Figure la is used for a right propagating wave and the stencil in Figure lb is 
used for a left propagating wave. The upwind bias nature of these stencils is thus clearly evident. 
References [3], [2], [6], [7], [8] clearly show that the LBS is second-order accurate. 


Note that the last two terms in (14) and (15) represent the electric and magnetic loss (or source) 
terms. A key element in developing an accurate LBS scheme is proper treatment of these source 
terms. Several methods for treating these source terms have already been investigated [6], [7], 
[8] for aeroacoustic applications. For electromagnetics, several methods for treating these source 
terms were investigated during this project. For example, one can index both source terms in (14) 
and (15) at time level n, which Thomas [6] has shown to be unstable. Another example would be 
to index both source terms at time level n + 1 resulting in a semi-implicit method. Other examples 
involve applying an exponential transformation to (14) and (15) to eliminate the source term and 
then perform a linearization of the exponential terms in the discretized equations. However, the 
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Figure 1: One-dimensional upwind leapfrog computational stencils for right-going (a) and left- 
going (b) characteristics. 


method found to be most efficient and accurate is to index the self source term in (14) (i.e. P) at 
time level n + 1 and to index the coupled source term Q at time level n. This avoids a matrix 
solution at each grid point, and the formulation easily limits to the perfect conductor condition as 

(T oo. 


Using the stencils shown in Figure 1 and the source term indexing scheme outlined above, the 
resulting finite difference equations for (14) and (15) are 
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These equations can be rewritten in the form 
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where v — cAt/Ax is the Courant number. We now rewrite equations (18) and (19) as 
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Equations (20) and (20) are the update equations for lossy dielectric and magnetic materials. Note 
that as <r oc, then we have the PEC condition that P" +1 = Q" +l = 0 as required. 
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3 Heterogeneous Dielectrics 

One of the difficulties with the conventional FDTD algorithm is the error in treatment of mate- 
rial discontinuities. Recent research efforts have attempted to reduce this error source by suitable 
averaging of material properties across the interface or by interpolation or extrapolation of the 
electromagnetic fields near these material boundaries [11], [12]. The advantage of the LBS is that 
the characteristic based nature of the algorithm leads to a very natural treatment of dieletric inter- 
faces. Since the LBS works with characteristic variables, the slope of characteristic curves in each 
material will be different, and the physical boundary conditions permit an elegant and efficient 
implementation of a dielectric interface boundary condition. This numerical boundary condition 
implements the physics exactly, with no averaging, interpolation or extrapolation required. 

To implement the dielectric material interface boundary condition, consider the one-dimensional 
grid shown in Figure 2. The dielectric interface is located at grid point i and the dielectric materials 



Figure 2: One-dimensional computational grid for the LBS showing characteristic variables, a di- 
electric interface located at cell i, and corresponding field components and characteristic variables 
used for the surface boundary condition. 

can be lossy. The characteristic variables at grid point i. Pi and Q ir are split into two components 
each: P \ , Q \ , P 2 and Q 2 . The terms Pi and Q\ exist just to the left of the material interface as shown 
in Figure 2. The remaining terms P 2 and Q 2 exist just to the right of the material interface. For 
material 1, equation (20) is used to predict the value for P" +1 at the boundary and for material 
2, equation (21) is used to predict the value for Q" +l . To complete the implementation, the Q” +l 
and P 2 +1 terms must be updated. These terms are updated by enforcing the physical boundary 
conditions on the electromagnetic field at the material boundary. We can then solve for Q1 +i and 
P 2 n+I in terms of the "known" characteristic variables Pf* +l and Q 2 +l . To develop this procedure, 
the electromagnetic boundary conditions on the tangential field components are given by 

Byi = E f2 = (24) 
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For the right-going wave, substituting (24) and (25) into (8) gives 



p» +1 = if;*' + - 1 - p;' • 


(26) 

= ^(/ r , +« r , ) + ^('? +1 

-< 35 +l ) 

(27) 

Similarly, substituting (24) and (25) into (9) yeilds 
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Since P” +I and Q" +l are determined at boundary point i from the usual update equations (we 
treat them as "known" variables), it is necessary to express P 2 " f 1 and Q \ l+ 1 in terms of these 
variables. Rearranging (27) and (29) gives 
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where T ,. 2 and T \ ■> are reflection and transmission coefficients given by 
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From (30), it is clear that the right-going wave in material 2 is a sum of a transmitted portion 
of the right-going wave in material 1 plus a reflected portion of the left-going wave in material 
2. A similar argument can be made for the left-going wave in material 1. In fact, the reflection 
coefficients T , 2 can be shown to be identical to the classical Fresnel reflection coefficients. The 
transmission coefficients also have the same form as the Fresnel transmission coefficients. Special 
care needs to be taken when the LBS calculates the solution at grid points i - 1 and i + 1 for a 
material interface at grid point i. At grid point i - 1, the term Q? +l in (23) becomes Q\ l . At grid 
point i, the terms P/ 1 and Q\ l in (22) and (23) become P, n and Q%, respectively. At grid point i + 1, 
the term P’ , _ l in (22) becomes Pj'. Rearranging equations (18) and (19) for grid point i we have 
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where u\ = c\ At/ Ax, and // 2 = c>At / A.r. The terms ax, <i->, bx, b 2 refer to the a and b coefficients in 
(12) and (13) for materials 1 and 2, respectively. These equations are now easily solved for P," + 1 
and Q " + 1 and then (30) and (31) are applied to obtain P]' +l and Q 2 +x ■ 
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4 Outer Boundary Condition 

The outer radiation boundary condition is used to terminate the computational lattice and per- 
mit outgoing waves to pass unreflected through the lattice boundaries [13]. The FDTD algorithm 
uses a spatial central difference operator where it uses field values from neighboring cells to up- 
date solution variables. Thus it cannot be used at the terminating faces of the problem domain. 
For example, the solution for a wave propagating left to right will eventually require a grid point 
outside the domain. To terminate the computational lattice, an additional equation (boundary 
condition) is needed to solve the system and this introduces information into the solution that is 
not required by Maxell's equations. 

On the contrary, the LBS requires no extraneous boundary condition. For the present LBS im- 
plementation, like the Method of Characteristics [9], the interior point algorithm calculates the 
left-going characteristic at the left boundary (i.e. i = 0) and the right-going characteristic at the 
right boundary (i.e. i = irnax). Thus for the LBS, at grid point i = 0, equation (21) calculates C( 0 ) 
and the incoming right-going characteristic, /'(()), is specified as a boundary condition. This same 
analysis applies at the right boundary where (20) calculates P{imax) and the incoming left-going 
characteristic, Q(irnax), is specified as a boundary condition. Shang [14] has noted for charac- 
teristic based multidimensional and nonuniform grid problems, in principle, the local coordinate 
system can be rotated to align with the characteristics, and the compatibility equations provide an 
exact boundary condition. However, a simple, yet effective approximation for multidimensional 
characteristic based approaches is to set the incoming flux or characteristic variables at the outer 
boundaries to zero and let the interior point algorithm predict the outgoing variables. When the 
wave motion is aligned with a coordinate axis, this boundary condition is exact. 


5 Fourier Analysis 


Various excellent Fourier analyses of the LBS have already been completed [3], [2], [6], [7], 
[8]; therefore, only the important results and conclusions from these previous analyses will be 
reviewed in this report. Most of the information presented is summarized from [2]. The stability 
condition for the LBS is v < 1, where v is the Courant number v = cAt/Ax. The normalized phase 
velocity for the classical leapfrog or FDTD algorithm is given by 


c* z/sin(A:Ax) 
c sin ) 


(38) 


Making the substitution kAx = 2i x/N and uAt = 2nv/N, where N is the number of points per 
wavelength, (38) becomes 
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The leading error term of the phase speed error is given by 
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error in phase speed 


It is clear from this error term that the FDTD algorithm is second-order accurate. For the LBS, the 
dispersion relation is 

c? _ (2i/~ l)sin(7r/JV) 
c sin (2m; /N — ■n/N) 

The leading error term of the phase speed error for the LBS is 
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Results are shown in Figure 3 for the normalized phase speeds at different Courant numbers (v) for 
the FDTD method and the LBS. To achieve less than 1% phase speed error requires about N = 15 




Grid resolution ( cells /wavelength) Grid resolution (cells/wavelength) 


Figure 3: Percentage error in phase speed versus grid resolution for (a) the FDTD method; and (b) 
the LBS. Plot parameter is v, the Courant number. 

for the FDTD method and about N = 6 for the LBS. Note that the LBS has zero dispersion error 
for u = 1 and v - 0.5. Based upon these results, the LBS method is about 2-3 times as enonomical 
as the FDTD method for the same level of accuracy. 

6 Results 

Since the LBS has previously been applied to free-space propagation problems on uniform 
grids [4], [7], this report concentrates on one-dimensional model problems involving free-space 
propagation on nonuniform grids and reflection from lossy dielectric materials on both uniform 
and nonuniform grids. The problem space size is 1000 cells with nonperiodic boundary condi- 
tions. For the uniform grid, a space step size of 1 cm is used, the time step is 0.33 ps and the 
Courant number v = 1. For the boundary conditions, a Gaussian point source at i = 0 is used to 
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specify P( 0) and Q(1000) = 0. For many complex geometries, it is often desirable to implement 
nonuniform grids to reduce the computational effort and memory resources and to improve mod- 
eling accuracy. We define a nonuniform grid by using a mesh stretch ratio of M = Ax max / A x 1llin 
which is periodic every 10 cells. Figure 4 shows an expanded view of a typical one-dimensional 
non-uniform grid with a mesh stretch ratio of 3 and a periodicity of 10 cells. The Courant num- 


0 0.2 0.4 0.6 0.8 1 

x (meters) 


Figure 4: Section of a one-dimensional non-uniform grid with a mesh stretch ratio of 3 and a base 
cell size of 1 cm. The grid variation is periodic every 10 cells. 


ber for a nonuniform grid is defined by cAt/Ax m i n , where Ax rnu , is the smallest cell size in the 
nonuniform grid. 

The first problem is a free space propagation problem on a nonuniform grid with a mesh 
stretch ratio of 2 that is periodic every 10 cells. In this case, periodic boundary conditions are used 
and the Gaussian pulse is allowed to propagate for 724 meters, which leads to a time integration of 
90,504 time steps. The Courant number v = 0.8, the time step was At = 2.67 ns and the Gaussian 
pulse had a FWHM pulse width of 2.26 ns. This pulse contained significant spectral content up 
to 1 GHz. Figure 5 shows the error in the electric field after n = 90, 504 time steps for both the 
FDTD method and the LBS. Note for this particular problem, the error for the LBS is exceptionally 

T3 

rH 

CD 

•rH 
4-1 

U 
-H 

4-> 

V 

a) 

t~ H 

<u 

a 

•H 

u 

o 

u _ 
u 
a) 

cX° “ ^ • -* 

724726728730732734736738740 
x (meters) 



Figure 5: Percent error in electric field for a free space propagation problem on a nonuniform grid 
using the FDTD method and the LBS. 
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low. From further experimentation, it was demonstrated that the LBS provided excellent results 
(within 0.1% accuracy) up to a mesh stretch ratio of 3. 


The next problem involved reflection and transmission for a lossy dielectric half-space. A uni- 
form grid was used first with the dielectric half-space for 5 < x < 10 m with material parameters 
e r 2 = 4; — 0.02, <1 ^ = 0 and /x r 2 = 1 and with u = 1. This problem tests the implementa- 

tion of the dielectric surface boundary condition which exists in the grid at cell i - 500. Figure 
6 shows the time-domain scattered electric field for both the LBS and the FDTD method at cell 
i = 400. Note the agreement between the two methods is almost indistinguishable. The complex 



Time (ns) 


Figure 6: Time-domain electric field recorded at cell 1 = 400 for reflection from a lossy dielectric 
half-space using FDTD and the LBS on a uniform grid. 

reflection coefficient was calculated at cell i = 400 for both methods and is plotted in Figure 7 
along with the exact solution. The agreement between methods is again excellent. Next, a lossless 
dielectric was inserted in the uniform grid domain for 5 < x < 10 m with material parameters 
< v> = 80, a 2 = cr* 2 = 0 and /x r 2 = 1. The complex reflection coefficient was computed from the 
time-domain fields and is shown in Figure 8. Clearly the LBS is superior in this instance, with a re- 
flection coefficient that overlays the exact solution. Note that the LBS solution is exact for this case 
because of the exact implementation of the physical boundary conditions on the electromagnetic 
field for the LBS. The lossy dielectric results were not as accurate because the dielectric bound- 
ary conditions assumed a frequency-independent impedance. If the frequency dependent surface 
impedance were used, then the results would be more accurate. Finally, the same lossy dielectric 
problem was analyzed on a nonuniform grid which is periodic every 10 cells and with a mesh 
stretch ratio of 2. The material parameters were again e r2 = 4, a 2 = 0.02,^ = 0 and /z r2 = 1. 
The dielectric surface boundary condition remains at i = 500, which puts the x coordinate for 
the boundary at 7.24 meters. Electric field data was recorded at x = 5 meters (i.e. i = 346). The 
time-domain fields are shown in Figure 9 and Figure 10 shows the reflection coefficient magnitude 
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Figure 7: Reflection coefficient magnitude versus frequency for reflection from a lossy dielectric 
half-space using FDTD and the LBS on a uniform grid. 
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Figure 8: Reflection coefficient magnitude versus frequency for reflection from a lossless dielectric 
half-space using FDTD and the LBS on a uniform grid. 
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results. We clearly see that the LBS is superior on a nonuniform grid, with a reflection coefficient 
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Figure 9: Time-domain electric field recorded at cell i = 346 for reflection from a lossy dielectric 
half-space using FDTD and the LBS on a non-uniform grid. 

accuracy level within 2%. Further results for uniform and nonuniform grids and perfect conduc- 
tors can be found in [15]. 

7 Conclusions 

This report has extended the Linear Bicharacteristic Scheme for computational electromagnet- 
ics to model homogeneous and heterogeneous lossy dielectric and magnetic materials and perfect 
conductors. It was demonstrated that the LBS has several distinct advantages over conventional 
FDTD algorithms. First, the LBS is a second-order accurate algorithm which is about 2-3 times as 
economical. The LBS can also be made to have zero dispersion error in certain instances. Second, 
the LBS provides a more natural and flexible way to implement surface boundary conditions and 
outer radiation boundary conditions by using characteristics and an upwind bias technique pop- 
ular in fluid dynamics. Third, the LBS provides more accurate results on nonuniform grids. The 
upwind biasing provides a more flexible generalization to unstructured grids. A dielectric surface 
boundary condition was implemented and results were provided for several one-dimensional 
model problems involving lossy dielectric materials and free space. The results indicate that the 
LBS is a superior algorithm for treatment of dielectric materials, especially its performance on 
nonuniform grids. Based upon these results, the LBS is a very promising alternative to a conven- 
tional FDTD algorithm for many applications. Extensions to two and three-dimensional problems 
should be straightforward. 
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x (meters) 

Figure 10: Percent error in reflection coefficient magnitude versus frequency for reflection from a 
lossy dielectric half-space using FDTD and the LBS on a non-uniform grid. 
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